home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / calcpam.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  3.9 KB  |  202 lines  |  [TEXT/R*ch]
  1. /*
  2.  * CALCPAM.C - Calculate a log-odds matrix for given PAM distance.
  3.  *
  4.  * Copyright David T. Jones, June 1992
  5.  *
  6.  * The basic command syntax is
  7.  *
  8.  *     calcpam datfile pamdist
  9.  *
  10.  * where pamdist is an integer >= 1
  11.  *
  12.  * References:
  13.  *
  14.  * Jones D.T., Taylor W.R. and Thornton J.M. (1992) "The rapid generation
  15.  * of mutation data matrices from protein sequences". Comput. Appl. Biosci.
  16.  * Vol. 8 No. 3, 275-282.
  17.  *
  18.  * Dayhoff, M. O., Schwartz, R. M. & Orcutt, B. C. (1978).
  19.  * A model of evolutionary change in proteins. In Atlas of Protein Sequence
  20.  * and Structure, Vol. 5, Suppl. 3, Ed. M. O. Dayhoff, pp. 345-352.
  21.  * Natl. Biomed. Res. Found., Washington.
  22.  *
  23.  */
  24.  
  25. /* Uncomment next line for PC TURBO C */
  26. /* extern unsigned _stklen = 30000; */
  27.  
  28. #include <stdio.h>
  29. #include <math.h>
  30.  
  31. char buf[160];
  32.  
  33. /*    One letter amino acid code    */
  34.  
  35. char *rescodes = "ARNDCQEGHILKMFPSTWYVBZX";
  36.  
  37. /* Normalized frequencies of occurrence */
  38.  
  39. double freq[20];
  40.  
  41. /* Relative mutabilities - not actually needed here. */
  42.  
  43. double mutab[20];
  44.  
  45. /* Exchange probabilities */
  46.  
  47. double    pam[20][20], relodds[20][20];
  48.  
  49. int npams;
  50.  
  51. #define roundint(x) ((x >= 0) ? (int)(x+0.5) : (int)(x-0.5))
  52. #define skiprems(ifp) while(fgets(buf, 160, ifp) && fgetc(ifp) != ' ');
  53.  
  54. /* Calculate Relatedness Odds Matrix */
  55. void 
  56. calcrel()
  57. {
  58.     int             i, j;
  59.  
  60.     for (j = 0; j < 20; j++)
  61.     for (i = 0; i < 20; i++)
  62.         relodds[i][j] = pam[i][j] / freq[i];
  63. }
  64.  
  65. /* Raise 1 PAM matrix to given integral power (>= 1). One day I'll
  66. speed this routine up by combining powers of 2. */
  67. void
  68. pampow(mat, power)
  69. double mat[20][20];
  70. int power;
  71. {
  72.     int             i, j, p;
  73.     double          mul[20][20], result[20][20], sum;
  74.  
  75.     for (i = 0; i < 20; i++)
  76.     for (j = 0; j < 20; j++)
  77.         mul[i][j] = mat[i][j];
  78.  
  79.     while (--power)
  80.     {
  81.         for (i = 0; i < 20; i++)
  82.             for (j = 0; j < 20; j++)
  83.         {
  84.         for (sum = p = 0; p < 20; p++)
  85.             sum += mul[i][p] * mat[p][j];
  86.         result[i][j] = sum;
  87.         }
  88.  
  89.     for (j = 0; j < 20; j++)
  90.         for (i = 0; i < 20; i++)
  91.         mat[i][j] = result[i][j];
  92.     }
  93. }
  94.  
  95. FILE *chkopen(name)
  96. char *name;
  97. {
  98.     FILE *f;
  99.  
  100.     f = fopen(name, "w");
  101.     if (!f)
  102.     {
  103.         fprintf(stderr, "Unable to open %s for writing!\n", name);
  104.         exit(1);
  105.     }
  106.     return f;
  107. }
  108.  
  109. void results()
  110. {
  111.     int i, j;
  112.     char fname[40];
  113.     FILE *ofp;
  114.  
  115.     sprintf(fname, "mp_%d.mat", npams);
  116.     ofp = chkopen(fname);
  117.     for (i = 0; i < 20; i++)
  118.     {
  119.     for (j = 0; j < 20; j++)
  120.         fprintf(ofp, "%8.5f ", pam[i][j]);
  121.     fprintf(ofp, "\n");
  122.     }
  123.     fclose(ofp);
  124.  
  125.     sprintf(fname, "ro_%d.mat", npams);
  126.     ofp = chkopen(fname);
  127.     for (i = 0; i < 20; i++)
  128.     {
  129.     for (j = 0; j < 20; j++)
  130.         fprintf(ofp, "%9.5f", relodds[i][j]);
  131.     fprintf(ofp, "\n");
  132.     }
  133.     fclose(ofp);
  134.  
  135.     sprintf(fname, "lo_%d.mat", npams);
  136.     ofp = chkopen(fname);
  137.     for (i = 0; i < 20; i++)
  138.     {
  139.     for (j = 0; j < 20; j++)
  140.         fprintf(ofp, "%6.1f", 10.0 * log10(relodds[i][j]));
  141.     fprintf(ofp, "\n");
  142.     }
  143.     fclose(ofp);
  144.  
  145.     sprintf(fname, "md_%d.mat", npams);
  146.     ofp = chkopen(fname);
  147.     for (i = 0; i < 20; i++)
  148.     {
  149.     for (j = 0; j < 20; j++)
  150.         fprintf(ofp, "%4d", roundint(10.0 * log10(relodds[i][j])));
  151.     fprintf(ofp, "\n");
  152.     }
  153.     fclose(ofp);
  154. }
  155.  
  156. main(argc, argv)
  157.     int             argc;
  158.     char           *argv[];
  159. {
  160.     int i, j;
  161.     FILE *ifp;
  162.  
  163.     if (argc != 3)
  164.     {
  165.     fprintf(stderr, "usage: calcpam datfile pamdist\n");
  166.     exit(1);
  167.     }
  168.  
  169.     sscanf(argv[2], "%d", &npams);
  170.     if (npams < 1)
  171.     {
  172.         fprintf(stderr, "PAM distance must be > 0!\n");
  173.         exit(1);
  174.     }
  175.  
  176.     ifp = fopen(argv[1], "r");
  177.     if (!ifp)
  178.     {
  179.         fprintf(stderr, "Unable to open datfile!\n");
  180.         exit(1);
  181.     }
  182.  
  183.     skiprems(ifp);
  184.     for (i=0; i<20; i++)
  185.         fscanf(ifp, "%lf", &freq[i]);
  186.  
  187.     skiprems(ifp);
  188.     for (i=0; i<20; i++)
  189.         fscanf(ifp, "%lf", &mutab[i]);
  190.  
  191.     skiprems(ifp);
  192.     for (i=0; i<20; i++)
  193.         for (j=0; j<20; j++)
  194.         fscanf(ifp, "%lf", &pam[i][j]);
  195.     fclose(ifp);
  196.  
  197.     if (npams > 1)
  198.     pampow(pam, npams);
  199.     calcrel();
  200.     results();
  201. }
  202.